package org.javanum.number;

import static org.javanum.number.Complex.COMPLEX_INFINITY;
/**
 * Copyright 2010 Scott Fines
 * <p/>
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * <p/>
 * http://www.apache.org/licenses/LICENSE-2.0
 * <p/>
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
 * either express or implied. See the License for the specific
 * language governing permissions and limitations under the License.
 * <p/>
 *
 * //TODO -sf- document!
 * Author: Scott Fines
 * Date: Mar 5, 2010
 */
public class ComplexMath {

    private ComplexMath(){}
    private static final Complex TWO_I = Complex.valueOf(0,2);
    private static final Real TWO = Real.valueOf(2);
    private static final PolarComplex POLAR_TWO_I
                                    = PolarComplex.valueOf(2,Math.PI/2);

    /*Sin functions*/

    /**
     * <p>Computes the sin of the passed in Complex number</p>
     *
     * 
     * @param z the number to compute the sin of
     * @return the sin of z
     */
    public static Complex sin(Complex z){
        if(Complex.isInfinite(z))
            throw new ArithmeticException(
                    "sin(infinity) is undefined.");
        if(Complex.isNaN(z))
            return Complex.NaN;

        Complex iZ = Complex.I.multiply(z);
        Complex iZNeg = Complex.I.negate().multiply(z);

        return exp(iZ).subtract(exp(iZNeg)).divide(TWO_I);
    }

    public static PolarComplex sin(PolarComplex z){
        if(PolarComplex.isInfinite(z))
           throw new ArithmeticException(
                    "sin(infinity) is undefined.");
        if(PolarComplex.isNaN(z))
            return PolarComplex.NaN;

        PolarComplex iZ = PolarComplex.I.multiply(z);
        PolarComplex iZNeg = PolarComplex.I.negate().multiply(z);

        return exp(iZ).subtract(exp(iZNeg)).divide(POLAR_TWO_I);
    }

/*----------------------------------------------------*/
    /*Arcsin functions*/
    
    public static Complex arcsin(Complex z){
    	/*
    	 * if z=a+ib, then
    	 * a must be restricted to the domain (-pi,pi),
    	 * where arcsin is valid. Otherwise, this
    	 * function is multivalued
    	 */
    	//TODO don't forget to restrict the domain!
    	Complex root = sqrt(Complex.ONE.subtract(pow(z,2)));
    	Complex diff = Complex.I.multiply(z);
    	Complex ln = ln(diff.add(root));
    	
    	return Complex.valueOf(ln.imaginary,-ln.real);
    }

    public static PolarComplex arcsin(PolarComplex z){
    	/*
    	 * TODO Don't forget to restrict the domain!
    	 */
    	PolarComplex root = sqrt(
    				PolarComplex.ONE.subtract(pow(z,2)));
    	PolarComplex diff = PolarComplex.I.multiply(z);
    	PolarComplex ln = ln(diff.add(root));
    	return ln.divide(PolarComplex.I);
    }

/*----------------------------------------------------*/
    /*Csc functions*/

    public static Complex csc(Complex z){
        return Complex.ONE.divide(sin(z));
    }

    public static PolarComplex csc(PolarComplex z){
    	return PolarComplex.ONE.divide(sin(z));
    }

/*----------------------------------------------------*/
    /*Cos functions*/

    public static Complex cos(Complex z){
        if(COMPLEX_INFINITY.equals(z))
            throw new ArithmeticException(
                    "cos(infinity) is undefined");
        Complex iZ = Complex.I.multiply(z);
        Complex iZNeg = Complex.valueOf(0,-1).multiply(z);

        return exp(iZ).add(exp(iZNeg)).divide(TWO);
    }

    public static PolarComplex cos(PolarComplex z){
        if(PolarComplex.isInfinite(z))
           throw new ArithmeticException(
                    "sin(infinity) is undefined.");
        if(PolarComplex.isNaN(z))
            return PolarComplex.NaN;

        PolarComplex iZ = PolarComplex.I.multiply(z);
        PolarComplex iZNeg = PolarComplex.I.negate().multiply(z);

        return exp(iZ).add(exp(iZNeg)).divide(TWO);
    }

/*----------------------------------------------------*/
    /*Arccos functions*/

    public static Complex arccos(Complex z){
        //TODO -sf- don't forget to bound the inputs!
    	Complex t = sqrt(pow(z,2).subtract(Complex.ONE));
    	
    	Complex r = ln(z.add(t));
    	
    	return Complex.valueOf(r.imaginary,-r.real);
    }

    public static PolarComplex arccos(PolarComplex z){
    	 //TODO -sf- don't forget to bound the inputs!
    	PolarComplex t = 
    		sqrt(pow(z,2).subtract(PolarComplex.ONE));
    	
    	PolarComplex r = ln(z.add(t));
    	return r.divide(PolarComplex.I);
    }

/*----------------------------------------------------*/
    /*sec functions*/

    public static Complex sec(Complex z){
        return Complex.ONE.divide(cos(z));
    }

    public static PolarComplex sec(PolarComplex z){
    	return PolarComplex.ONE.divide(cos(z));
    }

/*----------------------------------------------------*/
    /*tan functions*/

    public static Complex tan(Complex z){
        return sin(z).divide(cos(z));
    }

    public static PolarComplex tan(PolarComplex z){
    	return sin(z).divide(cos(z));
    }

/*----------------------------------------------------*/
    /*arctan functions*/

    public static Complex arctan(Complex z){
        //TODO -sf- bound inputs appropriately
    	Complex t1 = ln(Complex.I.subtract(z)).divide(TWO);
    	Complex t2 = ln(Complex.I.add(z)).divide(TWO);
    	
    	Complex r = t1.subtract(t2);
    	
    	return Complex.valueOf(r.imaginary,-r.real);
    }

    public static PolarComplex arctan(PolarComplex z){
        //TODO -sf- bound inputs appropriately
    	PolarComplex t1 = ln(PolarComplex.I.subtract(z)).divide(TWO);
    	PolarComplex t2 = ln(PolarComplex.I.add(z)).divide(TWO);
    	
    	PolarComplex r = t1.subtract(t2);
    	
    	return r.divide(PolarComplex.I);
    }

/*----------------------------------------------------*/
    /*cot functions*/

    public static Complex cot(Complex z){
    	return cos(z).divide(sin(z));
    }

    public static PolarComplex cot(PolarComplex z){
        return cos(z).divide(sin(z));
    }

/*----------------------------------------------------*/
    /*Exponential functions*/

    public static Complex exp(Complex z){
        if(Complex.isInfinite(z))
                return Complex.COMPLEX_INFINITY;
        if(Complex.isNaN(z))
                return Complex.NaN;

        double scale = Math.exp(z.real);
        double cos = Math.cos(z.imaginary);
        double sin = Math.sin(z.imaginary);

        return Complex.valueOf(scale*cos,scale*sin);
    }

    public static PolarComplex exp(PolarComplex z){
    	if(PolarComplex.isInfinite(z))
    		return PolarComplex.COMPLEX_INFINITY;
    	if(PolarComplex.isNaN(z))
    		return PolarComplex.NaN;
    	
    	double newMag = 
    			Math.exp(z.magnitude*Math.cos(z.angle));
    	double newAngle = z.magnitude*Math.sin(z.angle);
    	
    	return PolarComplex.valueOf(newMag,newAngle);    }

/*----------------------------------------------------*/
    /*logarithmic functions*/
    
    /**
     * Note that this method does NOT
     * satisfy ln(exp(z))==z, as z.imaginary may
     * be much larger than pi, while ln(z) will return
     * and imaginary part between -pi and pi. 
     * Therefore, the returned imaginary part of
     * ln(exp(z)) == z+2pi*n, for some integer n.
     * 
     *  //TODO -sf- finish documenting
     */
    public static Complex ln(Complex z){
    	double real = Math.log(z.modulus());
    	double imag = Math.atan2(z.imaginary,z.real);
    	
    	return Complex.valueOf(real,imag);
    }

    /**
     * Note that this method does NOT
     * satisfy ln(exp(z))==z, as z.imaginary may
     * be much larger than pi, while ln(z) will return
     * and imaginary part between -pi and pi. 
     * Therefore, the returned imaginary part of
     * ln(exp(z)) == z+2pi*n, for some integer n.
     * 
     *  //TODO -sf- finish documenting
     */
    public static PolarComplex ln(PolarComplex z){
    	if(PolarComplex.isInfinite(z))
    		return PolarComplex.COMPLEX_INFINITY;
    	if(PolarComplex.isNaN(z))
    		return PolarComplex.NaN;
    	if(PolarComplex.ZERO.equals(z))
    		throw new IllegalArgumentException
    				("Cannot compute ln(0)");
    	
    	if((Double.compare(Math.PI, z.angle)==0)||
    	   (Double.compare(-Math.PI, z.angle)==0))
    		throw new IllegalArgumentException
    			 	("Branch Cut encountered. " +
    			 	"Cannot compute the natural log" +
    			 	"of a negative real number");
    	
    	double x = Math.log(z.magnitude);
    	double y = z.angle<0?z.angle+(2*Math.PI):z.angle;
    	
    	return Complex.valueOf(x, y).polarForm();
    }

/*----------------------------------------------------*/
    /*Root functions*/

    public static Complex[] allRoots(Complex z, int root){
        //TODO -sf- implement!
        return null;
    }

    public static PolarComplex[] allRoots(PolarComplex z,int root){
    	PolarComplex[] roots = new PolarComplex[root];
    	double magnitude = Math.pow(z.magnitude, 1d/root);
    	double angleBase = z.angle/root;
    	double scale = Math.PI/root;
    	for(int k=0;k<root;k++){
    		double rootAdd = scale*k;
    		roots[k] = PolarComplex.valueOf(magnitude, angleBase+rootAdd);
    	}
    	return roots;
    }

    /**
     * Takes the first sqrt of z
     * @param z
     * @return
     */
    public static PolarComplex sqrt(PolarComplex z){
        //TODO -sf- deal with Complex.INFINITY
        double newMag = Math.sqrt(z.magnitude);
        double newAngle = z.angle/2;
        return PolarComplex.valueOf(newMag,newAngle);
    }

    public static Complex sqrt(Complex z){
    	double real = Math.sqrt(z.real+z.modulus())/Math.sqrt(2);
    	double imag = Math.sqrt(z.modulus()-z.real)/Math.sqrt(2);
    	imag*=Math.signum(z.imaginary);
    	return Complex.valueOf(real,imag);
    }

/*----------------------------------------------------*/
    /*Power functions*/
    public static Complex pow(Complex z, int pow){
    	return complexPow(z,Complex.valueOf(pow,0));
    }

    public static PolarComplex pow(PolarComplex z, int pow){
    	return 
    		PolarComplex.valueOf(Math.pow(z.magnitude,pow),
    							pow*z.angle);
    }

    public static Complex complexPow(Complex z, Complex w){
    	if(Complex.ZERO.equals(z)){
    		if(Complex.ZERO.equals(w)) 
    			throw new ArithmeticException("0^0 encountered");
    		else return Complex.ZERO;
    	}
    	
    	return exp(w.multiply(ln(z)));
    }
    
    public static PolarComplex 
    				complexPow(PolarComplex z, PolarComplex w){
    	if(PolarComplex.ZERO.equals(z)){
    		if(PolarComplex.ZERO.equals(w))
    			throw new ArithmeticException("0^0 encountered");
    		else return PolarComplex.ZERO;
    	}
    	return exp(w.multiply(ln(z)));
    }
}
